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Abstract 

The event-to-event fluctuations of hadron multiplicities are studied for a quark 
system undergoing second-order phase transition to hadrons. Emphasis is placed on 
the search for an observable signature that is realistic for heavy-ion collisions. It is 
suggested that in the 2-dimensional y-<p space the produced particles selected in a very 
narrow pr window may exhibit clustering patterns even when integrated over the entire 
emission time. Using the Ising model to simulate the critical phenomenon and taking 
into account a px distribution that depends on the emission time, we study in the 
framework of the void analysis proposed earlier and find scaling behavior. The scaling 
exponents turn out to be larger than the ones found before for pure configurations 
without mixing. The signature is robust in that it is insensitive to the precise scheme 
of simulating time evolution. Thus it should reveal whether or not the dense matter 
created in heavy-ion collisions is a quark-gluon plasma before hadronization. 



1 Introduction 

It is well known that a system undergoing a second-order phase transition exhibits large 
fluctuations and long-range correlations. If the symmetry of a quark system created in a 
heavy-ion collision is such that the quark-hadron phase transition is second order, then one 
would expect large fluctuations in the multiplicity of hadrons produced, not only from event 
to event (which happens normally in high-energy collisions), but also from one region to 
another in the geometrical space into which the particles are emitted at any given time. It is 
the latter type of fluctuations that we focus on as a possible signature of the phase transition, 
when an appropriate measure can be identified. 

The fluctuation of local hadron density implies the formation of spatial patterns that 
exhibit clusters of hadrons of varying sizes. There are also regions of no particles, which 
are referred to as voids. An investigation of the fluctuation of void sizes was carried out by 
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us recently in Ref. |L|. We found scaling behaviors, whose exponents were proposed as an 
observable signature of quark-hadron phase transition (PT). 

Recently an analysis of the NA49 data on Pb-Pb collisions at CERN-SPS was carried 
out according to the method described in Ref. [0. It is found || that the scaling behaviors 
exist under certain restricted conditions of the analysis, and that the scaling exponents are 
significantly lower than the ones predicted in [Ij. 

The void analysis studied in Jl| has its limitations. The most serious shortcoming is that 
it is an analysis of a system that is essentially static. The aim has been to find a measure 
suitable for quantifying effectively a given spatial pattern, or a set of them. However, the 
hadronization process in a heavy-ion collision is dynamical and takes place over an extended 
period of time. If the particles emitted throughout the whole duration are collected at every 
local region, then the patterns examined in the static scenario overlap one another, and the 
fluctuating patterns would disappear in the cumulative result. 

An analogy to this problem is the nature of rainfalls in two geographical regions with 
very different climates. It is not possible to judge from their annual precipitations (assuming 
that they are roughly the same) whether the rain in one region is stormy, or the other drizzly. 
It is necessary to make cuts in the time duration to learn about the spatial fluctuations in 
short time intervals. Similarly, in heavy-ion collisions it would be ideal if we could select a 
small At interval to analyze the spatial pattern of an event. That is, however, not feasible 
in a heavy-ion experiment. The only cut that can be made experimentally while keeping the 
y-(j) space for analysis is in the px variable, i.e., selecting particles in a small Apx interval. 
How to prepare the spatial configuration for our analysis in a way that corresponds to how 
the experimental data are to be analyzed is an essential part of the present work. 

As in which follows the theoretical framework or Ref. ||, we assume that the dense 
quark matter created in a heavy-ion collision is cylindrical in shape, expanding rapidly in the 
longitudinal direction and more gradually in the radial direction. It is further reasonable to 
assume that the thermal system is hot on the inside (with T > T c ) and cooler on the surface, 
where at T ~ T c the plasma undergoes a PT to hadrons. Thus the critical phenomenon 
occurs on a two-dimensional surface with the particle momenta being mainly normal to that 
surface. The fluctuation that we want to examine is in the hadron density from region to 
region in the y-(f) space for particles that are emitted in a small px interval. To simulate that 
fluctuation we make use of the theoretical result that the phase boundary between a first- 
order PT for some low quark masses and a crossover at higher quark masses is a second-order 
PT belonging to the same universality class as the Ising system [|], [5|. We map a section 
of the y-(f) space on the surface of the plasma cylinder created in a heavy-ion collision to a 
2D Ising lattice. We use the Ising model to simulate the spatial pattern of clusters of spins 
pointing in the same direction, and relate a cell with a net spin pointing in one direction 
to a hadron, and a cell with a net spin in the other direction to no-hadrons, i.e., locally 
and temporarily still in the quark state. That is the model that we have used in [3J, the 
essential part of which will be summarized in Sec. 2. 

Regarding what is simulated on the Ising lattice as a representation of the hadrons density 
fluctuation on the surface of the plasma cylinder at any given time, we have to relate that to 
the experimental situation of particles collected in a Apr cut, but accumulated throughout 
the entire hadronization time that can exceed 10 fm/c. It is a problem of preparing the 
temporally-integrated configuration on the lattice that can best approximate the real data. 
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The problem will be discussed in Sec. 3. The analysis of the mixed configurations that leads 
to scaling behaviors is the content of Sec. 4. 

2 Simulation and Analysis of Hadronic Density 
Fluctuation 

An essential feature of critical behavior is that the critical exponents of different systems 
belong to the same universality class if the systems have the same symmetry and dimension 
@. Our use of the Ising model to simulate the spatial patterns at PT is based on the 
recognition that for appropriate masses of the light and strange quarks the quark-hadron 
PT is in the same universality class as the Ising system |4|, |J. That is fortunate, since the 
simulation of Ising lattice is extremely simple. However, the realistic quark masses may 
not have exactly the values that correspond to the ones on the phase boundary where the 
second-order PT occurs. It is highly likely that the realistic point in the phase diagram is in 
the broad region of a crossover above the boundary line characterized by the Ising system. 
Fortunately again, a crossover can also be easily simulated on the Ising lattice by providing 
a preferred direction for the spins to be aligned. We summarize here the main points in the 
relationship between the Ising problem and the hadron density problem, and then outline 
the void analysis, the details of which can be found in 

For a 2D lattice of size L 2 with each site having spin Oj = ±1, we define the spin aligned 
along the overall magnetization m L = J2jeL 2 a j by 

Sj = sgn(m L )a j , (1) 

where sgn(mi) stands for the sign of mi,. We relate hadronization at location i to a cell of 
size e 2 with net spin c; positive, where 

Ci = J2 S 3 ' ( 2 ) 
3 

with the sum being over all sites in the cell. If c$ < 0, no hadronization takes place at i. 
Hadron density is defined by 

Pi = \<$e(c i ) , (3) 

where A is an unspecified factor relating the lattice spins to the number of particles per unit 
area in y-<f> space. The scaling behavior to be examined should, of course, be independent of 
the normalization factor A. 

Since lattice spins fluctuate from site to site when T is near T c , pi also fluctuates from 
cell to cell, either positive or zero. When averaged over all cells of a configuration, and then 
averaged over all uncorrelated configurations at any given T, (p) is a smooth function of 
T, descending rapidly near T c as T is increased [|1], |!|. It is a crossover, since the first few 
orders of derivatives exhibit no signs of discontinuities. The value of T c is determined by 
examining the scaling behavior of the normalized factorial moments F q J7|, and is found to 
be T c = 2.315 in the Ising model units of J/k B . 
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In practice we work with a lattice having L = 288 and cells of size e = 4. For scaling 
behavior we study the dependences of our measures on the bin size 5x5. The average hadron 
density in a bin is pb = v~ x Yh=i Pi, where v = (5/e) 2 . Near T c , pb fluctuates from bin to bin, 
especially for small 5. It is that fluctuation we want to characterize by quantitative measures. 
Our proposal is to represent p& in the 2D space by topographical maps with various discrete 
elevations po- F° r every given p we define a bin to be "empty" when 

Pb < Po ■ (4) 

and further define a void to be a collection of contiguous empty bins, connected by at least 
one side between neighboring bins. An analogy for this analysis is to flood a rough terrain by 
water and characterize the landscape profile by either land (nonempty bin) or water (empty 
bin). Clearly, by studying the areas of the submerged regions at different water levels, 
one can capture the topographical properties of the terrain. In our problem the maximum 
possible value of pi in a cell is (e 2 ) 2 = 256 in units of A. A range of p between 20 and 100 
have been used in our void analysis in |[J. 

For the void analysis let be the size of the Ath void in units of bins 

V k = J20(Po~p b ) , (5) 

where the sum is over all empty bins in the Ath void. Then, let x k = I4/M be the fraction 
of bins in the lattice occupied by V*., where M = (L/5) 2 is the total number of bins, and let 
m be the total number of voids in the configuration. We define the normalized G moments 
to be 

1 m I -| m 

(6) 

For every given configuration, G q can be determined for various bin sizes, i.e., as a function 
of M. G q can fluctuate widely from configuration to configuration. That fluctuation can be 
described by a probability distribution P(G q ). While many moments of that distribution 
can be studied, we have proposed the two lowest moments 

(G q )= JdG q G q P(Gq) , (7) 

S q = (G g \nG g ) , (8) 

the latter being the derivative with respect to p of the pth moment evaluated at p — 1 • 

For pure configurations simulated on the Ising lattice, scaling behaviors have been found 
for both (G q ) and 

(C,)ocM^ , (9) 
S q oc M a " , (10) 
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for a range of values of po 0- Furthermore, j q and a q depend on q linearly as 

l q = c + cq , 



(11) 



o q = s + sq . (12) 

Thus the values c and s (which are the slopes of slopes) are concise characterizations of the 
fluctuation behavior near the critical point. 

It should be stressed that although the 2D Ising model can be solved exactly for an infinite 
system[||, the conventional description of its critical behavior is in terms of the temperature 
T in the vicinity of T c . For a quark-hadron system created in a heavy-ion collision, T is 
not measurable, nor controllable. Our proposed measures of the critical behavior can be 
determined in high-energy nuclear experiments without any knowledge about T. 

3 Mixed Configurations for an Evolving System 

In the previous section we described the observable measures that can be applied to the 
heavy-ion data on the one hand, and can extract salient features of critical configurations on 
the other. Now we focus on the issue of preparing the configurations that correspond more 
closely to the realistic situation of a quark-gluon plasma expanding, cooling, and hadronizing. 

While a detailed modeling of the hydrodynamical expansion of a quark-gluon system is 
inappropriate for our purpose here, we can emphasize the main features of such an expanding 
system and incorporate them in our preparation procedure. Our assumption is that the 
temperature profile of the plasma cylinder is such that T is high (> T c ) in the interior and 
decreases monotonically toward the surface until T c , where quark-hadron PT takes place 
at the surface. The hadrons are emitted nonuniformly from the surface with transverse 
momenta pt that can be measured for each particle. The average (pt) depends on the radial 
pressure gradient, which in turn depends on time, since the longitudinal expansion lowers the 
temperature throughout the cylindrical interior and the radial expansion lowers the radial 
pressure gradient. Disregarding high-pr hard scattering, which affects (pt) only in a minor 
way, we may take the general hydrodynamical implication to be that (px) decreases with 
the pressure gradient, and therefore with the evolution time t. 

If the pt distribution at any given time were narrow, but when integrated over all times 
it gives the usual exponential damping in px for the event distribution dN/dpT, then our 
problem would be simple. In the idealized situation of a one-to-one correspondence between 
Pt and t, a narrow Apt cut in pt would result in a narrow At cut in emission time. That 
would provide us with a narrow temporal window to view the spatial pattern of hadronization 
activities on the cylinder surface, and the pure configurations simulated on the Ising lattice 
would suffice to model the effects of the critical behavior, as discussed in the previous section. 

The reality is, however, quite the opposite. At any given t the pt distribution is broad, 
so the particles emitted into a small Apt interval can come from a wide range of t. That 
introduces a complication related to the correlation between successive emission times. But, 
first, let us describe the pT-t correlation by a model formula that captures the essence of 
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the hydrodynamical character of the pr distribution at different times 0. We adopt the 
Gaussian formula 



P(pr,t) = N(t)exp{-[p T -p (t)] 2 /2pl(t)} 



(13) 



where 



p (t) = 2 ft GeV/c 
and N(t) is the normalization factor that preserves the condition 

f dp T P(p Tl t) = 1 



(14) 



(15) 



for all t > 1. It is clear that, as t increases, the Gaussian peak shifts toward lower pj- and 
its width also decreases correspondingly. If we make a narrow Ap T cut, say 0.3 < Pt < 0.35 
GeV/c, there are contributions from a wide range of t that are not negligible. 

We note that the accuracy of Eq. fljjD is not important for our purpose, since firstly 
the pt distribution depends on the evolution model and on many parameters in the initial 
condition, and secondly we use it only for simulating configurations at different times in the 
evolution process in the framework of the Ising model, which is a simple generator of PT 
configurations that ignore the realistic complications of surface irregularity, quark density 
fluctuations, surface tension, etc. However, Eq. (|13|) does contain the essence of overlapping 
Pt distributions at successive times, a feature that we want to incorporate in our simulation. 
The variable t in Eqs. ( fL3|) and ([14]) is used only in the sense of representing time in the 
simulation (to be detailed below), not as real time in heavy-ion collisions in some physical 
unit. We emphasize that the correspondences between the lattice variables and realistic 
quantities in heavy-ion experiments have never been fixed. In Eq. (^) we did not specify A 
that relates lattice spins to hadron density; the scaling behaviors we seek are independent 
of A. Here we similarly do not know how to specify the precise relationship between the 
number of steps in the evolution of the Ising lattice and the actual time of hydrodynamical 
evolution of a quark system. We now consider how many steps to take. 

If at any given t we have on the Ising lattice a configuration of hadrons and voids, then it 
is very reasonable to expect that at 1 fm/c later in real time the previous void-region will have 
a higher probability to hadronize than the region where hadronization has already occurred. 
Of course, there are fluctuations in local temperature, in surface curvature, pressure gradient, 
quark density, and other local properties of the statistical and collective system, all of which 
are unknown and can affect hadronization locally. Our procedure is primarily to let the 
Ising system evolve by updating the spin configuration using the Wolff logarithm ||10|| . Let 
S denote the number of steps of updating the Ising spin orientations that corresponds to 1 
fm/c in real time. Not knowing what S is, but realizing that the Ising configuration does 
not change drastically in one or two steps, we adopt the following two options. 
Option 1. S = 5 

After every five steps of updating the Ising configuration, we reverse the direction of what 
is "up" to allow the cells in the void region in the previous configuration to be preferred for 
hadronization. That is, denoting the number of Ising steps by nS, we let t — n in Eq. fllTf), 
and replace Eq. ([]]) by 




(16) 
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Note that if S were 1, and the Ising configuration changes very slowly, then Eq. ( ]T6[ ) applied 
to (H) and @ would suggest that p, L > and p, L < regions approximately alternate as n 
increases. That is an extreme scenario. By letting 5* = 5, we allow statistical readjustment 
of the lattice spins as time progresses. The (— l) n factor in (|16D is put in by hand to take 
into account the hydrodynamical feature that any quark region is likely to hadronize at a 
later time due to radial expansion and cooling. 
Option 2. S = 50 

After fifty steps it is likely that some of void regions will have evolved to contain Pi > 
cells. It is then not necessary to impose sign flip by hand, so we can continue to use Eq. 
(H) for Sj and rely on the fluctuations of the Ising spins to generate new configuration that 
correspond to successive emission time. That is, with the number of Ising steps being nS, 
we use t = n in Eq. (|13"D, and without the (— l) n factor in ([H]) we proceed with the evolution 
process. This is an option representing a scheme quite different from the one above. 

We now describe the algorithm for generating the temporally-integrated configurations 
(TIC) in either option of identifying t. At each t we simulate on the Ising lattice a pattern 
of hadrons and voids. Let us call that a frozen configuration. At the location of every cell 
where pi > 0, we assign a value of according to the distribution P(pr, t) given in Eq. fljjp. 
We choose to keep only those particles whose pr fall within the interval Ap? that we select, 
say 0.3 < pr < 0.35 GeV. The cell location of each of those selected particles are recorded 
in a new configuration, which is to become the TIC. Now let t increase by one step, repeat 
the above procedure, and enter more particles into the TIC at the cell positions where the 
selected hadrons are located. We continue to do this up to t m steps, when the accumulated 
particles at various places in the 2D space constitute our final TIC for that sum. The value 
of t m is determined by the requirement that the average hadron density for the whole lattice 
is approximately equal to that in the frozen case (no t evolution and any p?) at the same 
temperature T. In this way we generate a TIC, starting always with t = 2 so that po(t) 
ranges from 1 down to 0.1 GeV/c. 

For our analysis in the next section, we use a lattice of size L = 288, with cells of size 
e = 4, and bin sizes varying from 5 = 8 to 32. We select Ap? to be between 0.2 to 0.3 
GeV/c, and the corresponding t m turns out to be 9. We generate 5 x 10 3 configuration for 
each sample for analysis, using both options for t evolution. 

4 Scaling Behaviors of the Temporally-Integrated 
Configurations 

Having described how the TIC are prepared, we proceed directly to the analysis discussed 
in Sec. 2 and present the results here. 

Although the simulation can be done at any T in the vicinity of T c and the analysis can 
be done at any elevation po, we show only the scaling behaviors for T = T c and po = 40 [in 
units of A in Eq. (0)]. We shall consider both options of t evolution discussed in the previous 
section and label the corresponding figures by either S = 5 or S = 50. 

In Fig. 1(a) we show (G q ) vs M in log- log plot for S = 5. Clearly, it exhibits very good 
scaling behavior for all values of q considered. It means that voids of all sizes exist in the 
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TIC. Since the behavior is very similar to the one found in |1| for the pure or frozen con- 
figurations, the implication is that the mixing of configurations selected in the Apr window 
and integrated over the emission time does not lead to a homogenization of the clustering 
patterns. This behavior is possible only if each of the underlying (frozen) configuration ex- 
hibits critical fluctuation. Taking a small portion of each of such configurations at different 
t and adding them up retain the fluctuating nature of the void structure that gives rise to 
the scaling behavior. This is a very fortunate feature that makes possible the detection of 
the critical behavior in heavy-ion collision. 

From the slopes of the straight-line fits in Fig. 1(a), which are denoted by 7 9 in Eq. @, 
we show the dependence of 7 9 on q in Fig. 1(6). Evidently, Eq. (|TT|) is a good description of 
that linear dependence. The slope c of slopes is a numerical characterization of the scaling 
behavior of (G q ). Its value in this case is c = 0.92, which is significantly larger than the 
value of 0.7 found in |lj for pure (frozen) configurations. 

The same analysis is carried out in option 2 of t evolution. The corresponding results are 
shown in Figs. 2(a) and (b). The value of c = 0.97 is slightly larger, but is sufficiently close 
to the value from option 1 that one can reasonably infer the essential independence of the 
critical fluctuation on the procedure of simulating the temporal evolution. 

We have also investigated the dependences on p . In Fig. 3 we show the values of c for 
the two options in t evolution as functions of po- It should be understood that by giving the 
values of c we are asserting the existence of scaling such as those shown in Figs. 1(a) and 
2(a) for all the values of po examined. No essential difference between the two options of 
time evolution has been found. The feature in Fig. 3 to be noted is the wide range of po from 
40 to 100 in which the scaling behavior can be found. To give a sense of the normalization 
for those numbers, it should be remarked that the average (p) is around 78. That feature 
and the values of c ranging between 0.75 and 0.96 may now be regarded as the quantitative 
signature of quark-hadron phase transition that should be checked in heavy-ion collisions. 

For all that have been done for (G q ) above, we can do the same for S q defined in Eq. (§). 
Fig. 4(a) shows the scaling behavior of S q vs M, whose exponents o q are shown in Fig. 4(b), 
all for S — 5. The value of the slope s, defined in Eq. (|T2|) , is 0.82 and p = 40. For 5* = 50 
the corresponding figures are shown in Figs. 5(a) and (b), and the value of s is 0.89. Finally, 
the dependence of s on p is shown in Fig. 6. The range of po is the same as for (G q ), and 
the range of s is between 0.7 and 0.9. 

For a rough comparison we mention that in the real data of NA 49 the scaling behaviors 
of {G q ) and S q are found only in very restricted ranges of p and the values of c and s are 
significantly lower 0. 

For comparison, we have also performed the void analysis on randomly generated statis- 
tical configurations. We have chosen two methods to generate such configurations. 
Method 1 . We use random-number generator to assign a value of pi between and 160 at 
site i, for i ranging over all sites on the lattice. We choose that range because the average 
value of pi simulated in the foregoing at T c is around 78. We collect 5 x 10 3 such configurations 
and call them the random sample. 

Method 2. Generate first a configuration on the L x L sites of the lattice in the Ising 
model as in [1]. Choose randomly a cell i and transfer the hadron density pi there to a 
new configuration at the same cell location i. Repeat this procedure (L/e) 2 times, each 
time starting with a new uncorrelated configuration; furthermore, each time the transfer 
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of pj is made to the jth cell of the same new configuration. At the end of the (L/e) 2 th 
transfer, we have completed the construction of a new configuration, which is regarded as a 
mixed configuration. Then we start afresh and generate a new mixed configuration, adding 
to what we call the statistical sample. Specifically, we use L = 288 and e = 4. There 
are no correlations among the various cells in these configurations, although the underlying 
dynamics that generates the original configurations is that of the Ising model. 

The results of our void analysis on the random and statistical samples are very similar. 
Scaling behavior is not found until p is nearly as large as the mean density and then only 
for a narrow range of po between 60 and 80. This is very different from our results shown in 
Figs. 3 and 6 that have a wide range of po i n which the scaling behavior exists, beginning 
with po ~ (p) /2. 

Based on the above results we believe that dynamically generated critical behavior in the 
experimental data can be distinguished from the statistical background. 



5 Conclusion 

What we have demonstrated in this work is that critical fluctuation in quark-hadron phase 
transition can be observable in heavy-ion collisions, when appropriate measures are sought 
for in the data analysis. The details of how we have simulated the configurations are not 
as important as the insistence that a narrow Apx cut be made. There is no chance that 
the experimental data at RHIC can reveal any evidence of the pattern formation in hadron 
density if such a cut in p T is not applied. Scaling behavior is the canonical characteristics of 
critical phenomena, so to discover its existence in the data of heavy-ion collisions should be 
a primary area of investigation for the signatures of quark-gluon plasma. Since temperature 
is not directly measurable, and fluctuations in multiplicities do not manifest themselves 
naturally as in critical opalescence, our task is nontrivial and must rely on the identification 
of observables that are carefully constructed to reveal the scaling properties. Of course, if 
scaling is not found in the data, it could be that the phase transition is of first order, or 
is a crossover of such smoothness that no interesting fluctuations exist. The application of 
our method of analysis to the data can at least shed some light on the nature of the phase 
transition. 

There is a possibility that scaling behavior may exist even in the absence of a second- 
order PT of the type discussed here. Such a possibility has been proposed earlier [[II] as 



a manifestation of self-organized criticality \12\. Although the scaling discussed there is in 



terms of cluster sizes, there is no doubt that a void analysis can detect its signals, if it exists. 
The method proposed here should be sensitive to any hadronization process that does not 
produce particles uniformly in a short time interval. Application of the method to the RHIC 
data should be revealing whether or not scaling behavior can be found. 
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Figure Captions 

Fig. 1 Option 1 in time-evolution, S — 5: (a) Scaling behavior of (G q ), (b) Dependence 
of the scaling exponent 7 g on q. 

Fig. 2 Option 2 in time-evolution, S = 50: (a) Scaling behavior of (G q ), (b) Dependence 
of the scaling exponent 7 g on q. 

Fig. 3 Dependences of c on p for the two options. 

Fig. 4 Option 1 in time-evolution, S = 5: (a) Scaling behavior of S g , (b) Dependence of 
the scaling exponent a q on q. 

Fig. 5 Option 2 in time-evolution, 5 = 50: (a) Scaling behavior of S q , (b) Dependence 
of the scaling exponent a q on q. 

Fig. 6 Dependences of s on p for the two options. 
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